In [1]:
# a szokásos rutinok betöltése
%pylab inline
from scipy.integrate import * # az integráló rutinok betöltése
from ipywidgets import * # az interaktivitásért felelős csomag
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
from IPython.core.display import HTML
In [2]:
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit"
value="Click here to toggle on/off the raw code."></form>''')
Out[2]:
Felhasznált irodalom:¶
David C. Johnston: Thermodynamic Properties of the van der Waals Fluid, https://arxiv.org/abs/1402.1205
In [3]:
rc('text', usetex=True) # az abran a xticks, yticks fontjai LaTeX fontok lesznek
In [4]:
def nyomas(V,T):
pp = 8 *T/(3*V-1)-3/V**2
return (pp)
In [5]:
def Gibbs(V,T):
gV = -6/V + 8/3 * T/(3 * V-1)- 8/3* T * log(3*V-1)
return (gV)
In [6]:
def Freepot(V,T):
FV = -3/V - 8/3* T * log(3*V-1)
return (FV)
In [7]:
#Gibbs(.9,.2)
In [8]:
#Freepot(2,1)
A van der Waals-gáz izotermái, $p(V,T)$ görbe¶
In [9]:
Npont=200
T1=0.85
T2=0.9
T3=0.95
T4=1.
T5=1.05
T6 = 1.1
Vmin = 0.5
Vmax = 4
VV=linspace(Vmin,Vmax,Npont) #mintavételezési pontok legyártása
figure(figsize=(8,6))
plot(VV,nyomas(VV,T1),label=r'$T/T_c=$'+str(T1),lw=2,ls='-',color='brown')
plot(VV,nyomas(VV,T2),label=r'$T/T_c=$'+str(T2),lw=2,ls='-',color='blue')
plot(VV,nyomas(VV,T3),label=r'$T/T_c=$'+str(T3),lw=2,ls='-',color='green')
plot(VV,nyomas(VV,T4),label=r'$T/T_c=$'+str(T4),lw=2,ls='-',color='red')
plot(VV,nyomas(VV,T5),label=r'$T/T_c=$'+str(T5),lw=2,ls='--',color='black')
plot(VV,nyomas(VV,T6),label=r'$T/T_c=$'+str(T6),lw=2,ls='--',color='black')
legend(loc='upper right',fontsize=11)
title(r'$p(V)$ diagram',fontsize=20)
xlabel('V',fontsize=20)
ylabel('p',fontsize=20)
xlim(0.3,Vmax)
ylim(0,1.7);
xticks(fontsize=15)
yticks(fontsize=15);
grid();
In [10]:
def eqsVdW(p): # a ket egyensulyban levo fazis terfogatanak szamitasa
# a Gibbs-potencial es a nyomas azonos a ket fazis terfogatanal
V1, V2 = p
eq1=Gibbs(V1,T)-Gibbs(V2,T)
eq2=nyomas(V1,T)-nyomas(V2,T)
return (eq1, eq2)
Gibbs-potenciál, G = U-TS+pV van der Waals-gázra¶
In [11]:
Npont=200
T = 0.87
Vmin = 0.5
Vmax = 6
# A ket egyensulyban levo fazis terfogatanak szamitasa:
V1, V2 = fsolve(eqsVdW, (0.5, 2)) # a ket fazis terfogata
#print(V1,V2)
p0=nyomas(V1,T)
print ("A ket fazis terfogata V1 = ", V1, "és V2 = ",V2)
print ("check: ",eqsVdW((V1, V2)))
VV=linspace(Vmin,Vmax,Npont) #mintavételezési pontok legyártása
pp=nyomas(VV,T)
gk = Gibbs(Vmin,T)
gg=Gibbs(VV,T)
yy=gg-gk
#print(VV,res[0])
figure(figsize=(20,6))
subplot(1,2,1)
plot(VV,nyomas(VV,T),label=r'$T/T_c=\, $'+str(T),lw=3)
plot([V1,V2],[nyomas(V1,T),nyomas(V2,T)],lw=3,color="red")
# filled area from Maxwell construction
VV=linspace(V1,V2,Npont)
pfill=nyomas(VV,T)
#fill_between(VV, pfill, p0, where=pfill >= p0, facecolor='', hatch="\\")
fill_between(VV, pfill, p0, where=pfill >= p0, facecolor='', hatch='/')
fill_between(VV, pfill, p0, where=pfill <= p0, facecolor='', hatch='/')
legend(loc='upper right',fontsize=25)
title(r'$p(V)$ diagram',fontsize=30)
xlabel(r'$\hat{V}$',fontsize=30)
ylabel(r'$\hat{p}$',fontsize=30)
xticks(fontsize=25)
yticks(fontsize=25)
xlim(0.35,3)
grid()
ylim(0,1);
#savefig('/home/cserti/Dropbox/python/Termo/pV_VdW_abra_Maxwell.eps'); # Abra kimentese
subplot(1,2,2)
plot(pp,yy,lw=3)
title(r'$G(T,p)$ Gibbs-potenci\'al', fontsize=30)
xlabel('p',fontsize=30)
ylabel('G',fontsize=30);
xticks(fontsize=25)
yticks(fontsize=25);
xlim(0.15,1.)
ylim(-1.5,-0.4);
#print("T = ",T)
print("p_0 = ",p0)
Maxwell-konstrukció¶
In [12]:
Npont=200
T = 0.85
Vmin = 0.5
Vmax = 6
# A ket egyensulyban levo fazis terfogatanak szamitasa:
V1, V2 = fsolve(eqsVdW, (0.5, 2)) # a ket fazis terfogata
#print(V1,V2)
p0=nyomas(V1,T)
print("T = ",T)
print ("A ket fazis terfogata V1 = ", V1, "és V2 = ",V2)
print("p_0 = ",p0)
print ("check: ",eqsVdW((V1, V2)))
VV=linspace(Vmin,Vmax,Npont) #mintavételezési pontok legyártása
pp=nyomas(VV,T)
#print(VV,res[0])
figure(figsize=(8,6))
plot(VV,nyomas(VV,T),label=r'$T/T_c=\, $'+str(T),lw=3)
plot([V1,V2],[nyomas(V1,T),nyomas(V2,T)],lw=3,color="red")
# filled area from Maxwell construction
VV=linspace(V1,V2,Npont)
pfill=nyomas(VV,T)
#fill_between(VV, pfill, p0, where=pfill >= p0, facecolor='', hatch="\\")
fill_between(VV, pfill, p0, where=pfill >= p0, facecolor='', hatch='/')
fill_between(VV, pfill, p0, where=pfill <= p0, facecolor='', hatch='/')
legend(loc='upper right',fontsize=20)
title(r'$p(V)$ diagram',fontsize=20)
xlabel(r'$\hat{V}$',fontsize=20)
ylabel(r'$\hat{p}$',fontsize=20)
xticks(fontsize=15)
yticks(fontsize=15)
xlim(0.35,4)
grid()
ylim(0,1);
#savefig('/home/cserti/Dropbox/python/Termo/pV_VdW_abra_Maxwell.eps'); # Abra kimentese
Szabadenergia, F= U-TS potenciál van der Waals-gázra¶
In [13]:
Npont=200
T=0.85
Vmin= 0.35
Vmax = 11
# A ket egyensulyban levo fazis terfogatanak szamitasa:
V1, V2 = fsolve(eqsVdW, (0.5, 2)) # a ket fazis terfogata
#print(V1,V2)
p0=nyomas(V1,T)
print ("A ket fazis terfogata V1 = ", V1, "és V2 = ",V2)
print ("check: ",eqsVdW((V1, V2)))
Fk = Freepot(Vmin,T)
F1=Freepot(V1,T)-Fk
F2=Freepot(V2,T)-Fk
#print("V1 = ",V1,"F1 = ",F1)
#print("V2 = ",V2,"F2 = ",F2)
VV=linspace(Vmin,Vmax,Npont) #mintavételezési pontok legyártása
figure(figsize=(20,6))
subplot(1,2,1)
plot(VV,nyomas(VV,T),label=r'$T/T_c=\,$'+str(T),lw=3)
plot([V1,V2],[nyomas(V1,T),nyomas(V2,T)],lw=3,color="red")
# filled area from Maxwell construction
VV=linspace(V1,V2,Npont)
pfill=nyomas(VV,T)
#fill_between(VV, pfill, p0, where=pfill >= p0, facecolor='', hatch="\\")
fill_between(VV, pfill, p0, where=pfill >= p0, facecolor='', hatch='/')
fill_between(VV, pfill, p0, where=pfill <= p0, facecolor='', hatch='/')
legend(loc='upper right',fontsize=25)
title(r'$p(V)$ diagram',fontsize=30)
xlabel('V',fontsize=20)
ylabel('p',fontsize=20)
xticks(fontsize=25)
yticks(fontsize=25)
xlim(0.33,5)
grid()
ylim(0,1);
VV=linspace(Vmin,Vmax,Npont) #mintavételezési pontok legyártása
FF=Freepot(VV,T)
yy=FF-Fk
subplot(1,2,2)
plot(VV,yy,lw=3)
plot([V1,V2],[F1,F2],color="red", lw=3.0)
title(r'$F(T,V)$ szabadenergia',fontsize=30)
xlabel('V',fontsize=30)
ylabel('F',fontsize=30);
xticks(fontsize=25)
yticks(fontsize=25)
xlim(0,5);
grid()
ylim(-5,-2);
#print("T = ",T)
print("p_0 =",p0)
#tight_layout();
p-V diagram és a G Gibbs-potenciál¶
In [14]:
def rajzGibbs(Vmin,Vmax,V0,Npont,T,szin):
#mintavételezési pontok legyártása
V=linspace(Vmin,Vmax,Npont)
pp=nyomas(V,T)
gk = Gibbs(Vmin,T)
gg=Gibbs(V,T)
yy=gg-0*gk
figure(figsize=(20,6))
subplot(1,2,1)
plot(V,pp,lw=3)
p0=nyomas(V0,T)
plot(V0,p0,'o',markersize=15)
xlabel('V',fontsize=30)
ylabel('p',fontsize=30)
xticks(fontsize=25)
yticks(fontsize=25)
title(r'$p(V)$ diagram', fontsize=30)
xlim(0,Vmax)
grid(True)
subplot(1,2,2)
plot(pp,yy,color=szin,lw=3)
plot(p0,Gibbs(V0,T),'o',markersize=15)
xlabel('p',fontsize=30)
ylabel('G',fontsize=30)
xticks(fontsize=25)
yticks(fontsize=25)
title(r'$G(T,p)$ Gibbs-potenci\'al', fontsize=30)
grid(True) ;
In [15]:
def rajzFreepot(Vmin,Vmax,V0,Npont,T,szin):
#mintavételezési pontok legyártása
V=linspace(Vmin,Vmax,Npont)
pp=nyomas(V,T)
Fk = Freepot(Vmin,T)
FF=Freepot(V,T)
yy=FF-0*Fk
figure(figsize=(20,6))
subplot(1,2,1)
plot(V,pp,lw=3)
p0=nyomas(V0,T)
plot(V0,p0,'o',markersize=15)
xlabel('V',fontsize=30)
ylabel('p',fontsize=30)
title('$p(V)$ diagram', fontsize=30)
xticks(fontsize=25)
yticks(fontsize=25)
ylim(0,1.1)
grid(True)
subplot(1,2,2)
plot(V,yy,color=szin,lw=3)
plot(V0,Freepot(V0,T),'o',markersize=15)
xlabel('V',fontsize=30)
ylabel('F',fontsize=30)
xticks(fontsize=25)
yticks(fontsize=25)
title('$F(T,V)$ szabadenergia', fontsize=30)
ylim(-7,-4)
grid(True) ;
In [16]:
# Vmin = 0.57, Vmax = 5.3, T = 0.89, V0 --- csuszka
In [17]:
interact(rajzGibbs,
Vmin=FloatSlider(min=0.45,max=1,step=0.01,value=0.57,description='Vmin='),
Vmax=FloatSlider(min=2,max=10,step=0.1,value=5.,description='Vmax='),
V0 =FloatSlider(min=0.4,max=10,step=0.01,value=3,description='V0='),
T=FloatSlider(min=0.85,max=1.,step=0.01,value=0.89,description='T='),
Npont=IntSlider(min=200,max=200,step=10,description='Npoints='),
szin =Dropdown(options=['red','green','blue','darkcyan'],description='szín'));
p-V diagram és az F szabadenergia¶
In [18]:
# Vmin = 0.4, Vmax = 4.5, T = 0.86, V0 --- csuszka
In [19]:
interact(rajzFreepot,
Vmin=FloatSlider(min=0.35,max=1,step=0.01,value=0.4,description='Vmin='),
Vmax=FloatSlider(min=2,max=5,step=0.1,value=4.5,description='Vmax='),
V0 =FloatSlider(min=0.4,max=5,step=0.01,value=2,description='V0='),
T=FloatSlider(min=0.85,max=1.05,step=0.01,value=0.86,description='T='),
Npont=IntSlider(min=200,max=200,step=10,description='Npoints='),
szin =Dropdown(options=['red','green','blue','darkcyan'],description='szín'));
A van der Waals-gáz a kétfázisú tartománnyal¶
In [20]:
# coexistence line adatainak beolvasasa
# tablazatos formaban:
# T p Vf Vg
#
with open('/home/cserti/Dropbox/python/Termo/VanderWaals_coex_data.txt') as f:
sorok = f.readlines()
T_coex = []
p_coex = []
Vf_coex = []
Vg_coex = []
for sor in sorok[3:]:
T_coex.append( float(sor.split()[0]) )
p_coex.append( float(sor.split()[1]) )
Vf_coex.append( float(sor.split()[2]) )
Vg_coex.append( float(sor.split()[3]) )
In [21]:
Npont=200
Vmin= .5
Vmax = 4
T1=0.85
T2=0.9
T3=0.95
T3a= 0.975
T4=1.
T5=1.05
T6=1.1
tt=[T1,T2,T3,T3a,T4]
tfull = [T1,T2,T3,T3a,T4,T5,T6]
print("T = ",tfull) # ezeken a homersekleteken plottoljuk a p(V) gorbet
VV=linspace(Vmin,Vmax,Npont) #mintavételezési pontok legyártása
#figure(figsize=(8,6))
figure(figsize=(11,8))
figure(figsize=(8,6))
# p(V) izotermak a kritikus homerseklet folott
plot(VV,nyomas(VV,T5),label=r'$T/T_c=$'+str(T5),lw=2,ls='-',color='green')
plot(VV,nyomas(VV,T6),label=r'$T/T_c=$'+str(T6),lw=2,ls='-',color='green')
# a coexsistence line plottolasa
plot(Vf_coex,p_coex,color='blue',lw=2)
plot(Vg_coex,p_coex,color='blue',lw=2)
# T < Tc-re a p(V) plottolasa
for T in tt:
# A coexistence gorbenek az adatfile-ban megkeresni a T1,T2,T3,T3a,T4,T5,T6 homersekletu sorokat
T_indx=round((T-T1)/0.001)
# p(V)= allando a folyadek-goz fazis kozott
plot([Vf_coex[T_indx],Vg_coex[T_indx]],[p_coex[T_indx],p_coex[T_indx]],color='red',lw=2)
# p(V) a folyadek fazisra
VV=linspace(Vmin,Vf_coex[T_indx],Npont)
plot(VV,nyomas(VV,T),label=r'$T/T_c=$'+str(T),lw=2,color='black')
# p(V) a goz fazisra
VV=linspace(Vg_coex[T_indx],Vmax,Npont)
plot(VV,nyomas(VV,T),color='black',lw=2)
# p(V) instabil resz a folyadek-goz fazis kozott
VV=linspace(Vf_coex[T_indx],Vg_coex[T_indx],Npont)
plot(VV,nyomas(VV,T),color='black',ls='--',lw=2)
#legend(loc='upper right',fontsize=11)
title(r'$p(V)$ izoterm\'ak van der Waals-g\'azra',fontsize=20)
xlabel(r'$\hat{V}$',fontsize=20)
ylabel(r'$\hat{p}$',fontsize=20)
xlim(0.3,3.7)
ylim(0,1.5);
xticks(fontsize=15)
yticks(fontsize=15);
grid();
#savefig('/home/cserti/Dropbox/python/Termo/pV_VdW_abra_mod.eps'); # Abra kimentese
A folyadék-gőz fázist elválasztó p − T fázisgörbe van der Waals-állapotegyenletre¶
In [22]:
figure(figsize=(8,6))
title(r'$p(T)$ f\'azisdiagram',fontsize=20)
plot(T_coex,p_coex,'b-',lw=2)
plot([1],[1],'ro')
text(0.97,1.05,'Kritikus pont',fontsize=20,color='red')
text(0.83,0.75,r'folyad\'ek f\'azis',fontsize=20,color='black')
text(0.97,0.75,r'g\'az (g\H oz) f\'azis',fontsize=20,color='black')
xlabel(r'$\hat{T}$',fontsize=20)
ylabel(r'$\hat{p}$',fontsize=20)
xlim(0.8,1.1)
ylim(0,1.25);
xticks(fontsize=15)
yticks(fontsize=15);
grid();
#savefig('/home/cserti/Dropbox/python/Termo/pT_VdW_abra_mod.eps'); # Abra kimentese
A két fázis $V_g -V_f$ térfogatkülönbségének a $(T_c-T)/T_c$ dimenziótlan hőmérsékletől való függése¶
A közelítő függvény: $V_g-V_f=4 \sqrt{\epsilon} $, ahol $\epsilon = (T_c-T)/T_c$
In [23]:
figure(figsize=(8,6))
title(r'A k\'et f\'azis t\'erfogat\'anak k\"ul\"onbs\'ege, $V_g -V_f$',fontsize=20)
Tepsi=1-array(T_coex)
Vfg=array(Vg_coex)-array(Vf_coex)
plot(Tepsi,Vfg,'b-',lw=2)
# kozelito fuggveny a kritikus pont kornyeken
nn=90
plot(Tepsi[nn:],4*sqrt(Tepsi[nn:]), 'r-')
text(0.05,2.2,r'$V_g-V_f$',fontsize=22,color='blue')
text(0.05,2.0,r'numerikus ',fontsize=17,color='blue')
text(0.05,1.8,r'sz\'amol\'asb\'ol',fontsize=17,color='blue')
text(0.065,0.9,r'$V_g-V_f=4 \sqrt{\epsilon} $',fontsize=20,color='red')
text(0.065,0.7,r'k\"ozel\'{i}t\'es',fontsize=17,color='red')
xlabel(r'$\epsilon = (T_c-T)/T_c$',fontsize=20)
ylabel(r'$\hat{V}_g -\hat{V}_f$',fontsize=20)
xlim(0.,0.17)
ylim(0,3.);
xticks(fontsize=15)
yticks(fontsize=15);
grid();
#savefig('/home/cserti/Dropbox/python/Termo/Vt_VdW_abra_mod.eps'); # Abra kimentese
In [ ]: